# Cleaned jul17analysis.R - Complete code for all analyses in the article
# Set working directory and load libraries
library(tidyverse)
library(stargazer)
library(ggplot2)
library(dplyr)
library(ggpubr)
library(FSA)
library(patchwork)
library(reshape2)

# Load and prepare data
df <- read.csv2(file="DF_Jul17_4.csv")
names(df) <- tolower(names(df))

# Scale variables
df$ituc_gri <- df$ituc_gri - 7
df$ituc_gri <- df$ituc_gri * (-1)
df$ituc_gri <- scale(df$ituc_gri, center = FALSE)
df$yale_epi <- scale(df$yale_epi, center = FALSE)
df$wits_exports <- scale(df$wits_exports, center = FALSE)
df$wits_imports <- scale(df$wits_imports, center = FALSE)

# Filter data (no politicization data before 2004)
df <- subset(df, year >= "2004")
df$year <- df$year - 2003
df$year <- scale(df$year, center = FALSE)

# Create dataset groups
all <- df %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

countries <- subset(df, official_region=="none") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

no_countries <- subset(df, official_region !="none") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

# ========================================
# TABLE 2: MAIN REGRESSION MODELS
# ========================================

# Prepare models for regression (Table 1)
model_1 <- all %>% 
  rename(h1_norm_performance = yale_epi,
         h2a_market_exports = wits_exports,
         h2b_h2c_market_imports = wits_imports,
         h3_path_dependence = year,
         h4_h5a_politicization = env_pol,
         institutionalization = env_institutionalization)

model_2 <- no_countries %>% 
  rename(h1_norm_performance = yale_epi,
         h2a_market_exports = wits_exports,
         h2b_h2c_market_imports = wits_imports,
         h3_path_dependence = year,
         h4_h5a_politicization = env_pol,
         institutionalization = env_institutionalization)

model_3 <- all %>% 
  rename(h1_norm_performance = ituc_gri,
         h2a_market_exports = wits_exports,
         h2b_h2c_market_imports = wits_imports,
         h3_path_dependence = year,
         h4_h5a_politicization = lab_pol,
         institutionalization = lab_institutionalization)

model_4 <- no_countries %>% 
  rename(h1_norm_performance = ituc_gri,
         h2a_market_exports = wits_exports,
         h2b_h2c_market_imports = wits_imports,
         h3_path_dependence = year,
         h4_h5a_politicization = lab_pol,
         institutionalization = lab_institutionalization)

# Run main regressions (Table 1)
model_1_reg <- lm(institutionalization ~ h1_norm_performance + 
                    h2a_market_exports + 
                    h2b_h2c_market_imports +
                    h3_path_dependence +
                    h4_h5a_politicization, 
                  data = model_1)

model_2_reg <- lm(institutionalization ~ h1_norm_performance + 
                    h2a_market_exports + 
                    h2b_h2c_market_imports +
                    h3_path_dependence +
                    h4_h5a_politicization, 
                  data = model_2)

model_3_reg <- lm(institutionalization ~ h1_norm_performance + 
                    h2a_market_exports + 
                    h2b_h2c_market_imports +
                    h3_path_dependence +
                    h4_h5a_politicization, 
                  data = model_3)

model_4_reg <- lm(institutionalization ~ h1_norm_performance + 
                    h2a_market_exports + 
                    h2b_h2c_market_imports +
                    h3_path_dependence +
                    h4_h5a_politicization, 
                  data = model_4)

# Generate Table 1
stargazer(model_1_reg, model_2_reg, model_3_reg, model_4_reg,
          column.labels = c("env - all partners", "env - no countries", "lab - all partners", "lab - no countries"),
          omit = "Constant", 
          type="text", 
          out="July_19_multiple_regression.htm")

# ========================================
# TABLE 3: VARIANCE DECOMPOSITION (LMG METHOD)
# ========================================

# Install and load the relaimpo package for variance decomposition
if (!require(relaimpo)) {
  install.packages("relaimpo")
  library(relaimpo)
}

# Calculate relative importance using LMG method for all four models
lmg_1 <- calc.relimp(model_1_reg, type = "lmg", rela = TRUE)
lmg_2 <- calc.relimp(model_2_reg, type = "lmg", rela = TRUE)
lmg_3 <- calc.relimp(model_3_reg, type = "lmg", rela = TRUE)
lmg_4 <- calc.relimp(model_4_reg, type = "lmg", rela = TRUE)

# Extract the relative importance values
table3_data <- data.frame(
  Variable = c("h1_norm_performance", "h2a_market_exports", "h2b_h2c_market_imports", 
               "h3_path_dependence", "h4_h5a_politicisation"),
  Model_1 = round(lmg_1@lmg, 2),
  Model_2 = round(lmg_2@lmg, 2),
  Model_3 = round(lmg_3@lmg, 2),
  Model_4 = round(lmg_4@lmg, 2)
)

# Display Table 3
print("=== TABLE 3: VARIANCE DECOMPOSITION (LMG METHOD) ===")
print(table3_data)

# Create a formatted table using stargazer (optional)
stargazer(table3_data,
          type = "text",
          summary = FALSE,
          out = "Table_3_Variance_Decomposition.htm",
          title = "Table 3. Variance decomposition using Lindeman, Merenda, and Gold (1980) method",
          column.labels = c("Variable", "(1)", "(2)", "(3)", "(4)"),
          rownames = FALSE,
          digits = 2)

# Print summary information
cat("\n=== VARIANCE DECOMPOSITION SUMMARY ===\n")
cat("Model 1: Environmental - All Partners\n")
cat("Model 2: Environmental - No Countries\n") 
cat("Model 3: Labour - All Partners\n")
cat("Model 4: Labour - No Countries\n")
cat("\nMethod: Lindeman, Merenda, and Gold (1980) - LMG\n")
cat("Values represent relative importance (proportion of R²)\n")

# ========================================
# ROBUSTNESS CHECK I: WITHOUT POLITICIZATION
# ========================================

r_1 <- lm(institutionalization ~ h1_norm_performance + 
            h2a_market_exports + 
            h2b_h2c_market_imports +
            h3_path_dependence, 
          data = model_1)

r_2 <- lm(institutionalization ~ h1_norm_performance + 
            h2a_market_exports + 
            h2b_h2c_market_imports +
            h3_path_dependence, 
          data = model_2)

r_3 <- lm(institutionalization ~ h1_norm_performance + 
            h2a_market_exports + 
            h2b_h2c_market_imports +
            h3_path_dependence, 
          data = model_3)

r_4 <- lm(institutionalization ~ h1_norm_performance + 
            h2a_market_exports + 
            h2b_h2c_market_imports +
            h3_path_dependence, 
          data = model_4)

# Generate Robustness Check I table
stargazer(r_1, r_2, r_3, r_4,
          column.labels = c("env - all partners", "env - no countries", "lab - all partners", "lab - no countries"),
          omit = "Constant", 
          type="text", 
          out="July_19_multiple_regression_no_pol.htm")

# ========================================
# ROBUSTNESS CHECK II (SELF-CONTAINED)
# ========================================

# Create completely separate environment for Robustness Check II
{
  # Read fresh data for robustness check only
  df_rob2 <- read.csv2(file="DF_Jul17_4.csv")
  names(df_rob2) <- tolower(names(df_rob2))
  
  # Apply ONLY the ITUC transformation 
  df_rob2$ituc_gri <- df_rob2$ituc_gri - 7
  df_rob2$ituc_gri <- df_rob2$ituc_gri * (-1)
  
  # Filter for years >= 2004
  df_rob2 <- subset(df_rob2, year >= "2004")
  
  # Aggregate by trade agreement
  agg_df <- df_rob2 %>%
    group_by(trade.agreement) %>% 
    summarise(
      region_country = paste(unique(region_country), collapse = " "),
      yale_epi_sum      = mean(ifelse(is.na(yale_epi), 0, yale_epi)),  
      ituc_gri_sum      = mean(ifelse(is.na(ituc_gri), 0, ituc_gri)),  
      wits_imports_sum  = sum(ifelse(is.na(wits_imports), 0, wits_imports)),  
      wits_exports_sum  = sum(ifelse(is.na(wits_exports), 0, wits_exports)),
      avg_year          = median(ifelse(is.na(year), 0, year)),
      env_pol_avg       = median(ifelse(is.na(env_pol), 0, env_pol)),  
      lab_pol_avg       = median(ifelse(is.na(lab_pol), 0, lab_pol)),  
      env_legalization  = median(ifelse(is.na(env_legalization), 0, env_legalization)),
      lab_legalization  = median(ifelse(is.na(lab_legalization), 0, lab_legalization)),
      .groups = 'drop'
    ) %>%
    rename(
      env_performance = yale_epi_sum,
      labour_performance = ituc_gri_sum,
      imports = wits_imports_sum,
      env_politicisation = env_pol_avg,
      labour_politicisation = lab_pol_avg
    )
  
  # Print sample sizes
  cat("Total agreements:", nrow(agg_df), "\n")
  cat("Regional agreements (r only):", nrow(dplyr::filter(agg_df, region_country == "r")), "\n")
  
  # Run models
  rob2_env_1 <- lm(env_legalization ~ env_performance + imports + env_politicisation, data = agg_df)
  rob2_env_2 <- lm(env_legalization ~ env_performance + imports + env_politicisation, 
                   data = dplyr::filter(agg_df, region_country == "r"))
  rob2_lab_3 <- lm(lab_legalization ~ labour_performance + imports + labour_politicisation, data = agg_df)
  rob2_lab_4 <- lm(lab_legalization ~ labour_performance + imports + labour_politicisation, 
                   data = dplyr::filter(agg_df, region_country == "r"))
  
  # Generate table
  stargazer(rob2_env_1, rob2_env_2, rob2_lab_3, rob2_lab_4, 
            type = "text",
            out = "Robustness_Check_II_Clean.htm",
            omit = "Constant",
            title = "Robustness Check II",
            column.labels = c("-1", "-2", "-3", "-4"),
            dep.var.labels = c("Env_legalisation", "Lab_legalisation"),
            covariate.labels = c("Environmental Performance", "Labour Performance", "Imports", 
                                 "Environmental Politicisation", "Labour Politicisation"),
            keep.stat = c("n", "rsq", "adj.rsq", "ser", "f"),
            digits = 3,
            no.space = TRUE)
  
  # Print summaries
  print("=== ROBUSTNESS CHECK II RESULTS ===")
  print(paste("Total agreements:", nrow(agg_df)))
  print(paste("Regional agreements:", nrow(dplyr::filter(agg_df, region_country == "r"))))
  
  # Clean up robustness check variables
  rm(df_rob2, agg_df, rob2_env_1, rob2_env_2, rob2_lab_3, rob2_lab_4)
}

# ========================================
# REFRESH DATA FOR FIGURES 2 & 3  
# ========================================

# Reload and prepare main data to ensure clean state
df <- read.csv2(file="DF_Jul17_4.csv")
names(df) <- tolower(names(df))

# Reapply all transformations exactly as before
df$ituc_gri <- df$ituc_gri - 7
df$ituc_gri <- df$ituc_gri * (-1)
df$ituc_gri <- scale(df$ituc_gri, center = FALSE)
df$yale_epi <- scale(df$yale_epi, center = FALSE)
df$wits_exports <- scale(df$wits_exports, center = FALSE)
df$wits_imports <- scale(df$wits_imports, center = FALSE)

# Filter data
df <- subset(df, year >= "2004")
df$year <- df$year - 2003
df$year <- scale(df$year, center = FALSE)

# Recreate dataset groups
all <- df %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

countries <- subset(df, official_region=="none") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

no_countries <- subset(df, official_region !="none") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)



# ========================================
# CREATE REGIONAL SUBSETS FOR FIGURES 2 & 3
# ========================================

asean <- subset(df, official_region=="ASEAN") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

andean <- subset(df, official_region=="andean community") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

central <- subset(df, official_region=="central america") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

latin <- rbind(andean, central)

c_asia <- subset(df, official_region=="central asia") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

e_neighbor <- subset(df, official_region =="eastern neighborhood") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

s_neighbor <- subset(df, official_region=="southern neighborhood") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

balkans <- subset(df, official_region=="western balkans") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

w_africa <- subset(df, official_region=="west african states") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

sadc <- subset(df, official_region =="south african development community") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

eas <- subset(df, official_region=="eastern and southern africa") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

eac <- subset(df, official_region=="east african community") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

africa <- rbind(w_africa, sadc, eas, eac)

pacific <- subset(df, official_region =="pacific") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

carib <- subset(df, official_region=="caribbean") %>% 
  rename(env_institutionalization = env_legalization,
         lab_institutionalization = lab_legalization)

# Function to rename regions for figures
rename_regions <- function(data) {
  data %>%
    mutate(group = recode(group,
                          "countries" = "Countries",
                          "asean" = "ASEAN",
                          "latin" = "Latin America",
                          "c_asia" = "Central Asia",
                          "e_neighbor" = "Eastern Neighborhood",
                          "s_neighbor" = "Southern Neighborhood",
                          "africa" = "Africa",
                          "balkans" = "Western Balkans",
                          "carib" = "Caribbean",
                          "pacific" = "Pacific Island States"))
}

# ========================================
# PREPARE DATA FOR FIGURES 2 & 3 (FIXED)
# ========================================

# Function to rename regions for figures (using case_when instead of recode)
rename_regions <- function(data) {
  data %>%
    mutate(group = case_when(
      group == "countries" ~ "Countries",
      group == "asean" ~ "ASEAN",
      group == "latin" ~ "Latin America",
      group == "c_asia" ~ "Central Asia",
      group == "e_neighbor" ~ "Eastern Neighborhood",
      group == "s_neighbor" ~ "Southern Neighborhood",
      group == "africa" ~ "Africa",
      group == "balkans" ~ "Western Balkans",
      group == "carib" ~ "Caribbean",
      group == "pacific" ~ "Pacific Island States",
      TRUE ~ group
    ))
}

# Prepare data for Figures 2 & 3 (using dplyr::select)
env_data <- bind_rows(
  mutate(countries, group = "countries"),
  mutate(asean, group = "asean"),
  mutate(latin, group = "latin"),
  mutate(c_asia, group = "c_asia"),
  mutate(e_neighbor, group = "e_neighbor"),
  mutate(s_neighbor, group = "s_neighbor"),
  mutate(africa, group = "africa"),
  mutate(balkans, group = "balkans"),
  mutate(carib, group = "carib"),
  mutate(pacific, group = "pacific")
) %>%
  dplyr::select(value = env_institutionalization, group) %>%  # Use dplyr::select explicitly
  rename_regions()

lab_data <- bind_rows(
  mutate(countries, group = "countries"),
  mutate(asean, group = "asean"),
  mutate(latin, group = "latin"),
  mutate(c_asia, group = "c_asia"),
  mutate(e_neighbor, group = "e_neighbor"),
  mutate(s_neighbor, group = "s_neighbor"),
  mutate(africa, group = "africa"),
  mutate(balkans, group = "balkans"),
  mutate(carib, group = "carib"),
  mutate(pacific, group = "pacific")
) %>%
  dplyr::select(value = lab_institutionalization, group) %>%  # Use dplyr::select explicitly
  rename_regions()

# ========================================
# FIGURES 2 & 3: BOXPLOTS 
# ========================================

# Function to create Figures 2 & 3 (boxplots only)
create_boxplot <- function(data, variable_name) {
  # Kruskal-Wallis H Test
  kruskal_test <- kruskal.test(value ~ group, data = data)
  print(paste("Kruskal-Wallis test for", variable_name))
  print(kruskal_test)
  
  # Post-hoc Dunn's Test
  dunn_test <- dunnTest(value ~ group, data = data, method = "bonferroni")
  print(paste("Dunn's test for", variable_name))
  print(dunn_test)
  
  # Reorder the groups by descending means
  data <- data %>%
    group_by(group) %>%
    mutate(mean_value = mean(value, na.rm = TRUE)) %>%
    ungroup() %>%
    mutate(group = factor(group, levels = unique(group[order(-mean_value)])))
  
  # Boxplot to visualize the data
  boxplot <- ggboxplot(data, x = "group", y = "value",
                       color = "group", palette = "jco",
                       ylab = variable_name, xlab = "Regions") +
    scale_x_discrete(labels = function(x) {
      gsub(" ", "\n", x)
    }) +
    theme(axis.text.x = element_text(size = 10),
          axis.title = element_text(size = 12),
          legend.position = "none")
  
  # Save the plot
  plot_filename <- paste0(variable_name, "_boxplot.png")
  ggsave(plot_filename, plot = boxplot, width = 12, height = 8)
  
  return(boxplot)
}

# Generate Figures 2 & 3
env_plot <- create_boxplot(env_data, "Environmental_Institutionalization")
lab_plot <- create_boxplot(lab_data, "Labour_Institutionalization")

print(env_plot)
print(lab_plot)
# ========================================
# APPENDIX D: DETAILED HEATMAPS
# ========================================

# Load necessary libraries
library(ggplot2)
library(dplyr)
library(tidyr)
library(dunn.test)
library(reshape2)
library(patchwork)

# Load the data (using your exact approach)
df_filtered <- read.csv2("DF_Jul17_4.csv") %>% 
  filter(year >= 2004)

# Rename the columns (exactly as in your working code)
df_filtered_renamed <- df_filtered %>%
  rename(
    env_institutionalization = env_legalization,
    lab_institutionalization = lab_legalization
  )

# Filter out 'Central Asia' (exactly as in your working code)
df_filtered_renamed <- df_filtered_renamed %>%
  filter(region != "Central Asia")

# Perform Dunn's test for Environmental Institutionalization
dunn_env <- dunn.test(df_filtered_renamed$env_institutionalization, df_filtered_renamed$region, kw = TRUE, label = TRUE)

# Perform Dunn's test for Labour Institutionalization
dunn_lab <- dunn.test(df_filtered_renamed$lab_institutionalization, df_filtered_renamed$region, kw = TRUE, label = TRUE)

# Extract components from the list (ENVIRONMENTAL)
comparisons <- dunn_env$comparisons
p_adjusted <- dunn_env$P.adjusted

# Create a data frame
dunn_env_df <- data.frame(
  Comparison = comparisons,
  P.adjusted = p_adjusted,
  stringsAsFactors = FALSE
)

# Split the 'Comparison' column into 'Region1' and 'Region2'
dunn_env_df <- dunn_env_df %>%
  mutate(
    Region1 = sub(" - .*", "", Comparison),
    Region2 = sub(".* - ", "", Comparison)
  )

# Select relevant columns
dunn_env_df_selected <- subset(dunn_env_df, 
                               select=c("P.adjusted", "Region1", "Region2"))

# Pivot the data to wide format
dunn_env_df_pivoted <- dunn_env_df_selected %>%
  pivot_wider(names_from = Region2, values_from = P.adjusted, values_fill = list(P.adjusted = NA))

# Set 'Region1' as row names
dunn_env_df_pivoted <- tibble::column_to_rownames(dunn_env_df_pivoted, var = "Region1")

# Convert to matrix
dunn_env_matrix <- as.matrix(dunn_env_df_pivoted)

# Extract components from the list (LABOUR)
comparisons_lab <- dunn_lab$comparisons
p_adjusted_lab <- dunn_lab$P.adjusted

# Create a data frame
dunn_lab_df <- data.frame(
  Comparison = comparisons_lab,
  P.adjusted = p_adjusted_lab,
  stringsAsFactors = FALSE
)

# Split the 'Comparison' column into 'Region1' and 'Region2'
dunn_lab_df <- dunn_lab_df %>%
  mutate(
    Region1 = sub(" - .*", "", Comparison),
    Region2 = sub(".* - ", "", Comparison)
  )

# Select relevant columns
dunn_lab_df_selected <- subset(dunn_lab_df,
                               select=c("Region1", "Region2", "P.adjusted"))

# Pivot the data to wide format
dunn_lab_df_pivoted <- dunn_lab_df_selected %>%
  pivot_wider(names_from = Region2, values_from = P.adjusted, values_fill = list(P.adjusted = NA))

# Set 'Region1' as row names
dunn_lab_df_pivoted <- tibble::column_to_rownames(dunn_lab_df_pivoted, var = "Region1")

# Convert to matrix
dunn_lab_matrix <- as.matrix(dunn_lab_df_pivoted)

# Convert matrices to long format data frames and filter out NA values
dunn_env_long <- melt(dunn_env_matrix, varnames = c("Region1", "Region2"), value.name = "P.adjusted") %>%
  filter(!is.na(Region1) & !is.na(Region2))

dunn_lab_long <- melt(dunn_lab_matrix, varnames = c("Region1", "Region2"), value.name = "P.adjusted") %>%
  filter(!is.na(Region1) & !is.na(Region2))

# Define your custom ordering
y_axis_order <- c("Africa", "Caribbean", "ASEAN", "Countries", "Eastern Neighborhood", "Latin America")
x_axis_order <- c("Latin America", "Eastern Neighborhood", "Countries", "ASEAN", "Caribbean", "Western Balkans")

# Reorder the factor levels based on your custom order
dunn_env_long$Region1 <- factor(dunn_env_long$Region1, levels = y_axis_order)
dunn_env_long$Region2 <- factor(dunn_env_long$Region2, levels = x_axis_order)

dunn_lab_long$Region1 <- factor(dunn_lab_long$Region1, levels = y_axis_order)
dunn_lab_long$Region2 <- factor(dunn_lab_long$Region2, levels = x_axis_order)

# Define significance categories for environment
dunn_env_long <- dunn_env_long %>%
  mutate(Significance = cut(P.adjusted,
                            breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, 0.5, Inf),
                            labels = c("***", "**", "*", ".", "", "")))

# Define significance categories for labor
dunn_lab_long <- dunn_lab_long %>%
  mutate(Significance = cut(P.adjusted,
                            breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, 0.5, Inf),
                            labels = c("***", "**", "*", ".", "", "")))

# Custom color scale function
custom_fill_scale <- scale_fill_gradientn(
  colors = c("orange", "grey", "white"),
  values = scales::rescale(c(0, 0.1, 0.5, 1)),
  limits = c(0, 1),
  na.value = "white"
)

# Heatmap for Environmental Institutionalization with Significance and custom ordering
heatmap_env <- ggplot(dunn_env_long, aes(x = Region2, y = Region1, fill = P.adjusted)) +
  geom_tile() +
  geom_text(aes(label = Significance), color = "black") +
  custom_fill_scale +
  theme_minimal() +
  labs(title = "Environmental Institutionalization", fill = "P.adjusted") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "none")

# Heatmap for Labour Institutionalization with Significance and custom ordering
heatmap_lab <- ggplot(dunn_lab_long, aes(x = Region2, y = Region1, fill = P.adjusted)) +
  geom_tile() +
  geom_text(aes(label = Significance), color = "black") +
  custom_fill_scale +
  theme_minimal() +
  labs(title = "Labour Institutionalization", fill = "P.adjusted") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "none")

# Save and display heatmaps
ggsave("Appendix_D_Environmental_Heatmap.png", plot = heatmap_env, width = 10, height = 8)
ggsave("Appendix_D_Labour_Heatmap.png", plot = heatmap_lab, width = 10, height = 8)

# Display the heatmaps
print(heatmap_env)
print(heatmap_lab)

cat("\n=== APPENDIX D HEATMAPS COMPLETE ===\n")
cat("Custom ordering applied:\n")
cat("Y-axis: Africa, Caribbean, ASEAN, Countries, Eastern Neighborhood, Latin America\n")
cat("X-axis: Latin America, Eastern Neighborhood, Countries, ASEAN, Caribbean, Western Balkans\n")